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FOREWORD 


This is a final report on the research project, “Instability Wave Models and Supersonic Jet 
Noise.” The report titled “The Effects of Velocity Profiles on Supersonic Jet Noise” presented 
some of the earlier work done on this research project and, therefore, is not repeated here. Recent 
research activities were directed in the area of “Nonlinear Stability of Supersonic Jets.” 

The funding for this research was provided by the NASA Langley Research Center through 
the Grant NAG- 1-1 5 18. The grant was monitored by Dr. John M. Seiner of Aeroacoustics 
Branch (Fluid Mechanics and Acoustics Division), Mail Stop 165, NASA Langley Research 
Center, Hampton, VA 23681. 



NONLINEAR STABILITY OF SUPERSONIC JETS 


T. R. S. Shat* and S. N. TiwarP 
Department of Mechanical Engineering 
Old Dominion University, Norfolk, VA 23529-0247 

ABSTRACT 

The stability calculations made for a shock-free supersonic jet using the model 
based on parabolized stability equations are presented. In this analysis the large- 
scale structures, which play a dominant role in the mixing as well as the noise 
radiated, are modeled as instability waves. This model takes into consideration 
non-parallel flow effects and also nonlinear interaction of the instability waves. The 
stability calculations have been performed for different frequencies and mode num- 
bers over a range of jet operating temperatures. Comparisons are made, where ap- 
propriate, with the solutions to Rayleigh’s equation (linear, inviscid analysis with 
the assumption of parallel flow). The comparison of the solutions obtained using 
the two approaches show very good agreement. 
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1. INTRODUCTION 

The aerospace engineering community has placed a major emphasis on the 
development of a high speed civil transport aircraft. The success of this effort 
hinges on several factors. One of the major problems that needs to be addressed is 
the environmental issue of the noise generated by aircraft at take-off, cutback, and 
approach engine power settings. Jet noise represents a significant percentage of the 
overall noise generated. Hence, the reduction of jet noise is of utmost importance 
to future development of a supersonic civil aircraft. This has led to a revival of 
interest in developing techniques to predict and ultimately suppress noise radiated 
by supersonic jets. 

Several approaches have been developed for the prediction of jet noise. One 
such approach is the direct numerical simulation of the full Navier- Stokes equations 
for the estimation of the radiated noise. This technique, although exact, is not prac- 
tical because of the requirements of the computational resources. This has led to 
development of models based on simplified governing equations. Mankbadi et al. 1 
developed such a technique by using large-eddy simulation to model the flow field 
and calculated the far-field noise by applying Lighthill’s acoustic analogy. These 
techniques also require extensive computational power. Recently, there has been 
lot of focus on extending computational fluid dynamics techniques for aeroacoustic 
applications. The development of these techniques are not simple as the require- 
ments on numerical dissipation, dispersion, etc. for acoustic analysis are very rigid, 
see Tam 2 . 

A different approach for the prediction of jet noise is based on combining the- 
oretical and analytical methods. It is well established that turbulent flows possess 
coherent large-scale structures. These structures control the turbulent mixing pro- 
cess and also play a dominant role in the noise generation process. It has been 
successfully shown, in free shear layers such as plane mixing layers and jets, that 
the large-scale motions can be modeled as a random superposition of instability 
waves. The noise generated by these waves compare well with the experimental 
data. A good review of this approach is presented in Tam 3 . This model can be 
used to determine the two important characteristics of an instability wave, phase 
velocity and growth rate, which determine the mixing as well as the noise radiated. 

Most of the earlier theoretical studies of supersonic jets have considered linear, 
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inviscid analysis and have assumed the mean flow to be parallel. Tam and Burton 4 
took into consideration the non-parallel flow effects by using the method of multiple- 
scales. The real physical process of jets involve nonlinear interaction among all 
turbulence scales. The nonlinearities considered in this paper may be divided into 
nonlinear interaction between the instability waves and the mean flow, and the 
nonlinear self interaction of the instability waves. The former leads to a change in 
the mean-flow distribution, which in turn will modify the axial development of the 
instability waves. This kind of nonlinear effect can be taken into account by using 
measured mean velocity data. However, this data is not always readily available for 
different flow conditions. The effects of nonlinear self-interaction of the instability 
waves have been neglected in all the earlier studies by assuming that the amplitude 
of the disturbances, which initiate the instability waves near the nozzle exit, are 
small. This assumption may not be valid in all cases. 

In this study, nonlinear stability analysis of supersonic jets is considered. This 
analysis also accounts for non-parallel flows. The nonlinear interaction of the waves 
are predicted by the solving the parabolized stability equations (PSE). Bertoloti et 
al. 5 and Chang et al. 6 have applied PSE based technique for the stability analysis 
of incompressible and compressible boundary layers. The mean flow varies slowly 
in the axial direction enabling one to decompose the perturbation quantities into a 
rapidly varying wave-like part and a slowly varying shape function. The PSE can 
then be derived from the Navier-Stokes equations by applying parabolizing approx- 
imation to the shape function. The simplified governing equations can be solved 
by marching along the streamwise direction. This method is not computationally 
expensive. 

In this report stability calculations of supersonic jets using the PSE are pre- 
sented. The calculations have been performed for different frequencies and mode 
numbers over a range of jet temperatures, from cold to hot. In section 2, the de- 
velopment of the PSE model is presented and in section 3 computational procedure 
is described. The numerical results obtained for different cases are presented in 
section 4. Finally, in section 5, a summary of the results is given and also ways of 
improving the model and plans for future work are discussed. 



3 


2. ANALYSIS 

Consider the development of instability waves generated in a shock-free jet. It 
is assumed that the instability waves are governed by the compressible equations of 
motion, energy and state for a perfect gas. The flow variables are decomposed into 
a mean value and a perturbation quantity 

f(r,6,x,t ) = f(r,9,x ) + f'(r,0,x,t ) (2.1) 


where / is any flow variable and (r,9,x) is the polar coordinate system with x in 
the axial direction. The governing equations for the disturbances are obtained by 
substituting the form given by equation (2.1) in the governing equations and sub- 
tracting the equations for the mean flow. The resulting equations can be expressed 


4*4*4 * 4 *« 
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ox 1 ox or 
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where (f> is the disturbance vector and is given by <j) = (p 1 ,u' ,v' ,w' ,p') T . The 
coefficient matrices, T, A, B, Ve x and Voe are composed of a linear part and a 
nonlinear part. The linear part (denoted by superscripts /) contains only mean flow 
quantities and the nonlinear part (denoted by superscripts n) contains disturbance 
quantities , i.e. T = T* + T n , etc. Equation (2.2) can then be rearranged as 
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The governing equations (2.3) have to be parabolized to enable a marching 
solution in the streamwise direction. This is accomplished by decomposing the dis- 
turbances into a rapidly varying wave-like part and a slowly varying shape function. 
The wave part retains its ellipticity while the governing equations for the shape func- 
tion are parabolized. We assume that the disturbance vector <j> for an instability 
wave with a frequency mu and azimuthal mode number n can be expressed as 


m = oo n=oo 


<t>(r,6,x,t) = E E ^ mn (r, X^Xmn ($* 0 

m = — oo n= — co 

where the wave-like part Xmn is given by 

Xmn = expji a mn (t)d£ + n0 - mut j 


(2.5) 


( 2 . 6 ) 


ckmn is the axial wavenumber and ^ mn is the shape function vector given by ^ mn = 
(/3, ti, u, w , p ) T . As the shape function is assumed to vary slowly with the streamwise 
direction r, the second derivative d 2 'h/dx 2 is very small and is neglected. The 
governing equations for the shape function of a single mode (ra,n), using the form 
given by equation (2.5), can then be written as 


Dmn^mn + A * 


d<H, 


dx 


d* mn = yl d 2 * 


or 

+ F mn /A 


dr 2 


(2.7) 


where 


Amn — GXp 


. Jx 0 


The nonlinear forcing function F mn can be evaluated from the Fourier series expan- 
sion of F n 

m — oo n=oo 

E E A- exp [i(nd — mut)] (2.8) 


l= — OO 71— — OO 


In equation (2.7), the elliptic effect associated with the wave part is absorbed in 
matrices £), A and C and does not contribute to the upstream influence. However, 
the term dp/d: r allows upstream influence in subsonic regions of the flow. In order 
to make equation (2.7) truly parabolic, dp/dx is multipled by a constant Q given 

by 

f 1 Ml 

Q = { T+PFi JM?, 

1 


M x < 1 
M x > 1 


(2.9) 
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where M x is the local Mach number. This parabolizing procedure yield solutions for 
the parabolized Navier- Stokes equations which compare well with those obtained 
by the full Navier-Stokes equations, see Rubin 7 . 
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3. NUMERICAL SCHEME 


The parabolized stability equation (2.7) is solved numerically by finite- 
differencing the derivatives with respect to x and r. In the formulation given below, 
'Iq j represents the shape function with an axial location index i and radial location 
denoted by index j . The streamwise derivative at (i,j) is replaced by a second-order 
backward difference given by 

d<V 

= (3*.,, - Wi-u + Vi- 2 , j)/2Ax (3.1) 

for all axial locations except for the starting plane where a first-order backward 
differencing is used. The resulting discretized equation for the t-th axial plane is 

V itj =i,, J (4^ l _ lii -^_ 2ii )/2Ax + F/>t (3.2) 


D + ^-A + C-^-K^ 

2Ax dr dr 2 


where for convenience, the subscript mn which identifies the mode (m, n) is dropped. 
The radial derivatives are discretized using a fourth-order central differencing and 
axe given by 

^l\Ti \T/ . . . . 1 8lT/ . • . - RlLr • • . _L 11/ « 

(3.3) 


dV_ _ —^ 1,^+2 + 8Vjj+i ~ + ^»,. 7-2 

dr 


12Ar 

d 2 # -Vj,j+2 + 16^j,j-n ~ 30^j,j + - Vjj-2 4 s 

dr 2 ~ 12A r 2 K ' ’ 

The above finite differencing scheme is replaced by a second-order scheme for the 
grid point next to the boundary. The second-order scheme can be written as 


dV 

dr 


= (Vij+i -4,j + ,)/2Ar 


(3.5) 


d 2 V 

- 2*,,; + )/2Ar 2 (3.6) 

The finite difference formulations given above are substituted into equation (2.7). 
This coupled with the boundary conditions results in a block penta-diagonal system 
of equations (block size of 5 x 5) at each axial location for the unknown shape 
functions, V mn . 

The boundary conditions have to be prescribed for the eigenfunctions at the jet 
centerline and at a radial location of infinity. As the governing equations become 
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singular at r = 0, L’HospitaTs rule is applied at the centerline to obtain a new set of 
governing equations. The boundary conditions along the centerline takes a different 
form based on the azimuthal mode number. For example, for axisymmetric mode, 
i.e. n = 0, the boundary conditions at r =0 are given by 

u = v = 0, dxbjdr = dp/dr = 0 (3.7) 


Here, the velocity components u,v,w are in the r,0,x directions, respectively. The 
boundary condition on density is obtained from (3.7) and the discretized continuity 
equation. At the outer boundary, far from the centerline, all the shape functions 
Eire set equal to zero. For the nonlinear analysis, the far-field boundary conditions 
can be directly applied for all modes except for modes with zero frequency (m = 0) 
and axisymmetric mode (n = 0). The case of m = 0 and n = 0 is referred to as 
the mean flow correction mode. For this mode, the condition u = 0 is replaced by 
du/dr = 0. 

The solution procedure is as follows: start with a known solution for a mn and 
at x = xo, assume a value for a mn at the next axial location, march the PSE 
to the next station by solving for The new value of a mn is determined from 


d$ r 


&mn, new — ^ mn , old 




dx 


(3.8) 


and this process is repeated until the difference between the assumed and computed 
values of a mn is less than a given tolerance. The updating procedure given above is 
equivalent to normalizing the shape function such that d'b/dx is zero at a particular 
radial location. As the shape function vector \F is a function of r and contains five 
dependent variables, there is no unique wavenumber as in the liner theory for parallel 
flows. However, as shown by Chang et al. 6 , the total growth rate appears to be very 
weakly dependent upon the normalization. Further details of the normalization and 
its effect on the computed wavenumber are given in Chang et al. 6 . In all the results 
presented here, we have chosen to compute the wavenumber based on streamwise 
velocity using 

u*du/dx dr 

*’ JT I fi I 2 dr 


a 


mn, new 


— old 


(3.9) 


where * denotes complex conjugate. 
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4. NUMERICAL RESULTS 


The PSE has been used to analyze the linear stability of a supersonic jet is- 
suing from a circular nozzle operating at design Mach number of 2.0. It has been 
assumed that the axial evolution of the instability waves can be treated as an invis- 
cid phenomena and F mn in equation (2.7) is zero. The results obtained for parallel 
flows using the PSE are compared with the solutions to equation governing linear, 
inviscid instability wave, i.e., the compressible Rayleigh’s equation. The solutions 
to the Rayleigh’s equation (RES) are obtained using the method of Seiner, Bhat 
and Ponton 8 . For both the approaches, i.e., PSE and RES, the mean flow field is 
assumed to be known. It has been shown that the mean velocity profile can be 
closely approximated by a half- Gaussian function given by 


W(x) 


W c (x) r <h 

kE c (x)exp[— ln2(r — h)/b 2 ] r > h 


(4.1) 


where h(x ) is the radius of the potential core, b(x) is the half-width of the shear 
layer and W c (x) is the jet centerline velocity. For the present calculations, the 
values of these parameters axe obtained from computational data of Viswanathan 
et al. 9 . 

Calculations have been performed for a range of frequencies and the modes 
considered are n = 0, the axisymmetric mode, n = 1 and 2, the helical modes. 
These calculations were carried out over a range of jet operating temperatures. 
However, for conciseness, the results for n — 1 mode are shown here for a cold 
jet. The two most important characteristics of an instability wave are the phase 
velocity and growth rate. The phase velocity is defined as as 2n f j a r where / is the 
frequency and a r is the real part of the complex wavenumber a. The growth rate is 
defined as — cx t where a, is the imaginary part of the wavenumber. Figures 1 and 2 
compare the phase velocity and growth rates obtained using the two approaches at 
two different axial locations. It can be seen that over the entire range of Strouhal 
number (St = 2 frj/Wj) considered, the two solutions are in very good agreement. 
Similar agreements axe also obtained for n = 0 and 2. Figures 3 and 4 show 
the comparison as a function of downstream distance for axisymmetric mode (n 
= 0) at St = 0.54 and the helical mode (n = 1) at St = 0.49. Once again, the 
agreement is excellent for both the phase velocity and growth rate. As mentioned 
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before, calculations have also been performed for hot jets obtaining good agreement 
between the solutions using the two approaches. 

The results presented so far have been obtained for parallel flows. We now 
consider the solution of linear PSE with non-parallel flow effects. As in any marching 
scheme, an initial solution is needed to start the computations. To provide a good 
starting solution, the multiple-scales analysis can be used which accounts for non- 
parallel flow effects. It has been shown by Chang et al. 6 that using any incorrect 
initial conditions result in a transient phenomena close to the jet exit. Furthermore, 
it was observed that the PSE approach is capable of recovering the correct solution 
in the marching process provided the input initial conditions are close to the correct 
solution. In this study, the linear parallel stability equations are solved and the 
eigenvalue and eigenfunctions are used as the starting solutions. 

The linear PSE with non-parallel flow effects have been solved for several cases. 
Figures 5 and 6 show the results for two of these cases. Once again, calculations 
are carried out for a cold supersonic jet operating at its design Mach number of 
2.0. In figure 5, the variation of the wavenumber (real part) and the growth rate 
with axial distance is presented for axisymmetric mode (n = 0) at 5< = 0.54. 
The results from parallel stability theory are included for comparison. As can be 
seen, the wavenumber and the growth rates start off being less than those based 
on parallel theory. After some initial transients, both the wavenumber and growth 
rate are larger than those from parallel theory. This implies that the wave based 
on non-parallel flow has a lower phase velocity than the wave based on parallel 
flow. In addition, the wave corresponding to the non-parallel flow gets damped 
further downstream than that obtained using parallel flow assumption. The results 
for helical mode (n = 1) at St = 0.49 are presented in figure 6. Similar trends as 
in axisymmetric mode are observed. The only difference is that the wavenumber 
based on parallel and non-parallel mean flow are more or less the same. 

The eigenfunction distribution for the cases discussed above are presented next. 
Figure 7 shows the radial distributions of the eigenfunctions for the axisymmetric 
mode, St = 0.54 at X/Rj = 0. The eigenfunctions shown axe the radial velocity 
component, u, axial velocity component, u), and pressure, p. For this case, the az- 
imuthal velocity component, v, is zero. As can be seen, the axial velocity component 
is very large and peaks in the shear layer. These eigenfunctions further downstream, 
at X/Rj = 10, are presented in figure 8. The amplitude of the pressure eigenfunc- 
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tion has increased slightly while the amplitude of the axial velocity component has 
decreased substantially. The axial location for peak noise emission, based on the 
experimental observations of Troutt and McLaughlin 10 and Seiner, et al. 11 , occurs 
near each wave’s neutral point. This suggests that the amplitude of pressure would 
increase with downstream distance until the wave has become damped. 

The eigenfunctions for n = 1, Si = 0.49 at two axial locations, X/Rj = 0 and 
X/Rj = 10 are shown in figures 9 and 10, respectively. In this case, the azimuthal 
velocity component is nonzero. It can be seen that the trends are similar to that of 
the axisymmetric mode. The only difference is that the amplitude of the azimuthal 
velocity component, t), has increased. Converged solutions to the nonlinear PSE 
could not be obtained. Currently, this problem is being investigated. Some of the 
possible reasons for not getting converged solutions axe discussed in the next section. 
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5. DISCUSSIONS AND CONCLUSIONS 

In this study, a model for linear and nonlinear stability analysis of supersonic 
jets is developed using the parabolized stability equations. The linear PSE has 
been used to analyze the instability waves of a supersonic jet operating at design 
conditions, i.e., shock-free. For parallel flows, it has been shown that the PSE 
approach and the solutions to Rayleigh’s equation are in very good agreement over 
a range of frequencies and azimuthal mode numbers. The axial development of the 
instability waves also compared well. These calculations also considered a range of 
jet operating temperatures from cold to hot. The important characteristics of the 
instability waves, i.e., phase velocity and growth rate, are predicted well using the 

PSE. 

The linear PSE model has also been used to study the non-parallel flow effects. 
The computed solutions show that the growth rate is higher and the wave’s phase 
velocity is lower than that corresponding to the parallel flow. This expected trend 
is seen for instability waves of different frequencies and mode numbers. We have 
assumed that the dynamics of the axial development of the instability waves can be 
treated as an inviscid phenomena. At distances far downstream from the jet exit 
the instability wave would no longer grow in amplitude and is damped. At these 
locations the governing equations become singular and the integration contour must 
be deformed around the critical point, see Tam 12 and Tam and Morris 13 . In this 
study, the calculations are terminated close to the neutral point of the wave. 

The results presented here are preliminary and is by no means complete. In 
order to improve this model, there are several issues that need to be investigated. 
First of all, the solutions obtained from linear parallel stability equations have been 
used as the initial conditions for the marching scheme. The sensitivity of the initial 
conditions to the converged solutions have to be studied. Solutions obtained using 
the multiple-scales analysis, which accounts for non-parallel flow effects, can be 
used as the initial conditions for the eigenvalue and the eigenfunctions. Secondly, 
to circumvent the problem associated with the damped inviscid wave and to study 
the effects of Reynolds number, viscous terms have to be included in the analysis. 

As mentioned before, converged solutions to the nonlinear PSE could not be 
obtained. It appears that the initial amplitudes of the various waves considered 
play a critical role in getting converged solutions. Furthermore, the selection of the 
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waves, with nonzero initial amplitudes, also influences the iterative process. The 
effects of the initial amplitudes and the choice of the waves on getting converged 
solutions should also be investigated. The solutions obtained from nonlinear PSE 
analysis, after all the issues related to convergence have been identified and resolved, 
can be validated by comparing it to solutions computed from nonlinear numerical 
simulation, for example, the results of Viswanathan and Sankar 14 . 
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